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ABSTRACT 

We present a new efficient technique for measuring evolution of the galaxy luminosity 
function. The method reconstructs the evolution over the luminosity-redshift plane 
using any combination of three input dataset types: 1) number counts, 2) galaxy 
redshifts, 3) integrated background flux measurements. The evolution is reconstructed 
in adaptively sized regions of the plane according to the input data as determined by 
a Bayesian formalism. We demonstrate the performance of the method using a range 
of different synthetic input datasets. We also make predictions of the accuracy with 
which forthcoming surveys conducted with SCUBA2 and the Herschel Space Satellite 
will be able to measure evolution of the sub-millimetre luminosity function using the 
method. 
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1 INTRODUCTION 

A cornerstone of any working model of the formation of 
structure in the universe is knowledge of the galaxy lumi- 
nosity function (GLF). The GLF is a measure of the co- 
i moving space density of galaxies per interval in luminosity. 

■ Determining how the GLF changes with cosmological epoch 
, provides far-reaching insights into the processes that dictate 

■ the means by which galaxies form and evolve. This places 
' strong statistical constraints on evolutionary theories, en- 
abling determination of key characteristics such as the epoch 

' of galaxy formation, merger rates, transformations between 
population types and the history of the universe's global rate 
of formation of stars. 

The value of measuring the evolution of GLFs has been 
appreciated for several decades and as such, many tech- 
niques have been employed to achieve this aim. Without any 
loss of generality, the GLF can be expressed as the local GLF 
scaled by an evolution function. The simplest and most di- 
rect method of estimating this evolution function is a model 
independent one where the GLF is measured at different 
epochs and compared with the local GLF. This provides a 
set of discrete estimates of the evolution at specific epochs. 
However, adopting a model-based procedure yields signifi- 
cant advantages. A well-proven approach is to assume the 
evolution function depends only on luminosity and redshift 
(see next section) . For a given model of the evolution func- 
tion chosen a priori, the best fit model is determined, subject 
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to a set of observational constraints which contribute to dif- 
ferent regions in the luminosity-redshift {L — z) plane. The 
virtue of model-based techniques is that a model consistent 
with all available observations can be used to make predic- 
tions by extrapolation into regions where data are sparse or 
lacking. Furthermore, a model can be used to identify those 
datasets whose improvement would most efficiently increase 
the constraints on the evolution function. 

Model-based techniques are typically regarded as be- 
longing to one of two types. Parametric evolution models 
adhere to some preconceived notion regarding the evolution 
of galaxies, for example, pure lu minosity evolution or lumi- 
nosity -f density evolution (see IWalL Pope fc Scott 1 12008| . 
and references therein for recent examples). Although an ad- 
vantage of this type of modelling is that smooth functions 
are readily extrapolated, a major disadvantage is that the 
real evolution function may take on a very different form. 
So called 'free-form' methods are not biased in this way and 
allow for a greater degree of fiexibility over the L ~ z plane 
when attempting to determine the evolution function. 

Freefor m techniques date back over three decades. 
[Robertson ] 11978, 1980) developed an iterative freeform 
method to evaluate the evolution of the luminosity function 
for radio galaxies. The method was limited by allowing the 
evolution to vary only as a function of redshift (and not lumi- 
nosity) and not giving a satisfactory indica tion of the range 
of solu tions permitted by the observations. iPeacock fc Gull I 
(jTgU) introduced a new method that incorporated luminos- 
ity dependence and measured the uncertainty on the evolu- 
tion by considering the variation between different freeform 
model predictions. This improved technique saw successful 
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application (e.g., |Peacock|[l985l : iDunlop fc Peacock |[T99ol ') 
although by the authors' own admission, the full extent 
of the uncertainty still could not be verified by the range 
of models tested. Furthermore, the evolution function was 
modelled as a series expansion and the set of best-fit expan- 
sion coefficients determined by minimising using a non- 
linear search. This is both inefficient and does not guarantee 
that the global minimum has been found. 

In this paper, we propose a new method, inspired by 
Peacock & Gull but offering some significant improvements. 
The method computes the evolution as a discretised (i.e., 
pixellised) function over the L — z. plane. This has the ad- 
vantage that the evolution can be solved linearly, ensuring 
the best fit solution is always found with a single, highly 
efficient matrix inversion. Crucially, the full range of solu- 
tions permitted by the observations can be determined in a 
few simple extra computational steps. Despite being pixel- 
lised, we demonstrate how linear regularisation allows the 
evidence to be extrapolated into regions of the L — z plane 
that are lacking data or where data are sparse. Finally, we 
show how the pixellisation can be adapted according to the 
constraints provided by the data. 

The original motivation for this work was to investigate 
evolution of the luminosity function in the sub-millimetre 
(submm) waveband. The far-infrared (far-IR) and submm 
wavebands are particularly important for investigating the 
evolution of galaxies, not least because approximately half of 
the energy emitted by stars and active galactic nuclei since 
the big bang has been absorbed by dust an d then re-radiated 
in these wavebands (jPixsen et al. l[T998t l. In view of this, 
our interests lie in investigating not only the monochromatic 
luminosity functions at these wavelengths, but also the 'dust 
luminosity function', i.e., the space density of galaxies as a 
function of the total luminosity emitted by the dust in a 
galaxy. 

Presently, observational data in the submm are partic- 
ularly poor. However, this situation will begin to rapidly 
improve with the arrival of new submm instruments such 
as the second generation Submm Common User Bolometer 
Array (SCUBA2) and the Herschel Space Observatory {Her- 
scbel). Although the method we have developed is applica- 
ble to the luminosity function in any waveband, we describe 
and explore its use in this paper with an eye on its future 
application to forthcoming submm data. 

In section [2] we present the method. Section [3] demon- 
strates the method with a range of synthetic datasets provid- 
ing differing levels of constraints on the L — z plane. Finally, 
we summarise in section [4] and discuss practical aspects of 
applying the method to real data. 



2 THE METHOD 

The overall aim of the method is to provide an estimate of 
the evolution of the GLF. We therefore write the GLF, </!>, as 
the product of the local GLF, (jx) and an evolution function, 
E: 



<j,{L,z) = ML)E{L,z). 



(1) 



We have assumed that iJ is a bivariate function of total far- 
IR/submm luminosity, L, and redshift, z. We work in terms 
of the bolometric luminosity rather than a monochromatic 



luminosity, in keeping with the traditions adopted in submm 
astronomy. 

Until very recen tly, there wer e no direct submm mea- 
surements of (po (see lEales et al. 1 .2009'). Instead, the local 
GLF was estimated by extrapolating from shorter wave- 
length observations. In this way, the measurements had 
small statistical errors but, because of the uncertainty in the 
extrapolation, potentially large systematic errors. In what 
follows, we work towards determining E directly, accepting 
the fact that in application to real submm data, E poten- 
tially has systematic errors because of the uncertainty in 
(fio- Of course, it is a trivial exercise to re- write the formal- 
ism below in terms of (j}{L, z) instead of E[L, z) and use <j)o 
as an observational constraint. However, in order to deter- 
mine E from (j)(L, z) obtained in this way, one would have to 
renormalise by cj>o. This would therefore introduce the exact 
same systematic uncertainty on E as that encountered by 
computing E directly. 

The crux of the method is to write E as a, discrete func- 
tion, instead of the continuous form it has above. In this way, 
the evolution function takes on discrete values, Cj, in pixels 
which span the L — z plane. As we show below, this dis- 
cretisation turns the problem of determining E into a more 
desirable linear one. 

Observational constraints on E can be categorised into 
three different types: number counts, galaxy redshifts and 
integrated background flux. Number count data provide the 
number of detected galaxies per unit area of sky binned by 
flux. The number of galax;ies between two flux limits gives in- 
formation about the integral of the luminosity function over 
an extended region of the L — z plane, but does not provide 
information about whether the sources are luminous galax- 
ies at high redshift or galaxies with low luminosity at low 
redshift in that region. Writing this in terms of the lumi- 
nosity function, the expected number of detectable galaxies 
within a given flux bin i is 



dV{z) / dL(l>{L,z) 

z=0 Jl'Az) 



(2) 



where the lower and upper luminosity limits L] and 
correspond to the lower and upper flux limits of the bin 
and are therefore a function of redshift. These limits also 
depend on the frequency at which the flux is measured. 
Since L is a bolometric luminosity, a frequency depen- 
dent correction must be applied to obtain the monochro- 
matic rest-frame luminosity and hence observed flux 
Svivo) = (f 4- z)L,j{i') / A-Kd\, where v is the rest-frame fre- 
quency, Vo = v / {1 + z) IS the observer frame frequency and 
cLl is the luminosity distance. This correction is similar to 
a standard k-correction and depends on the galaxy spectral 
energy distribution (SED - see Section [3TT}- 

Using equations ((T| and ((2)1 , the number of galaxies that 
contribute to rii from a given pixel j in the L — z plane can 
be written 



AV{z) 



dL (j)o{L) — Cj aij 



(3) 



Here, the volume integral is written assuming the jih pixel 
in the L — z plane has fixed redshift boundaries. In general, 
the L — z plane can be divided up into pixels of arbitrary 
geometry in which case the integral would extend over the 
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portion of pixel j that contributes to bin i. The luminos- 
ity limits are capped by the luminosity range spanned by 
pixel j, giving rise to zero contribution if they lie outside 
the pixel. Note also that the luminosity limits have picked 
up an additional index, reflecting their dependence on the 
redshift limits of pixel j. The approximation in equation ^ 
arises due to the fact that the continuous evolution implicit 
in equation ^ has been replaced by the discrete evolution 
ej which is assumed constant over the entire pixel. The dis- 
cretised version of equation ((2| can thus be written: 



(4) 



where the sum acts over all pixels in the L — z plane although 
only a fraction will typically contribute to number count bin 
i. 

To incorporate redshift data, the procedure is essen- 
tially the same as for number counts. Redshift surveys pro- 
vide the most direct measurement of the luminosity function 
at a specific place in the L—z plane. Instead of binning galax- 
ies by flux, and evaluating the integral in equation ((2]) over 
^5 2 ^ oo, the data are binned by flux and redshift and the 
integral in equation ((2]) extends only over the redshift range 
of the bin. Instead of the quantities atj defined in equation 
(121, an analogous set of quantities bij are computed. 

Including integrated background fiux measurements re- 
quires a slightly different approach. Rather than compute 
the number of galaxies in a given data bin, the strength and 
spectral shape of the extragalactic background radiation in 
the far-IR/submm waveband puts a constraint on the inte- 
gral of the luminosity-weighted GLF over the entire L — z 
plane. The predicted integrated background flux, , as mea- 
sured at a given observer- frame frequency, Vo, is the sum of 
contributions from all galaxies over the whole L — z plane, 
i.e., 



dV{z) 



dL4>iL, z) s„{vo) ■ 



(5) 



As before, the flux Su{uo) is computed allowing for the SED 
dependent correction that converts bolometric luminosity to 
monochromatic luminosity in the rest frame of the galaxy. 
Writing this equation in its discretised form gives 

y"e, / dViz) / dLML)s^{uo,,) (6) 

j 

for the predicted integrated background fiux measurement 
at the ith observer-frame frequency Uo,i. 

From the observed galaxy number counts, n°, the ob- 
served number counts binned by redshift, n° i and the mea- 
sured background flux, (all of which are assumed to be 
statistically independent of each other), the statistic in 
terms of the discretised predicted quantities is 



^1 l„o 



^ -Y. ^1 + E ^ 

i = l i=JVi+l 
i = JVi+JV2 + l i 



In the above, there are A''i number count bins, number 
counts binned by redshift, Nt data points in total includ- 
ing the background flux measurements and the Oi are the 
1(7 measurement errors. To simplify, the last term deflnes 
the general quantity gi as an observed data point and pij 
refers to one of the corresponding quantities Oij, hij or dj. 
The set of values Cj that minimise are those that satisfy 
dx^ /dch = which gives 



9iPik/crf = BjPijPik/ai 



(8) 



By writing the vector element dk = Ei gi pik jal and the 
matrix element Muj = Tii Pij Pik / cr'i , the solution to the set 
of linear equations in ((Sjl can be written 



e = M"^d, 



(9) 



where e is a vector containing the elements Cj. 

In the presence of noise, the solution given by equation 
([9)1 is formally ill-conditioned and hence must be regularised. 
This is achieved by adding an extra term, the regularisation 
matrix R, weighted by the regularisation weight, A (see Sec- 
tion [23J: 



(10) 



e= (M-f- AR)"M. 



The corre sponding covariance matrix was derived by 
I Warren "fc Dye (2003) for this problem: 



C = N-AN(NR)^, N=(M + AR)"\ 



(11) 



This innocuous equation brings about the method's key ad- 
vantage over previous methods since the covariance matrix 
contains all of the information required to determine the un- 
certainty on the evolution function for a given set of model 
parameters (e.g.. A). The only extra step required to obtain 
the total error is to propagate the additional error that arises 
as a result of the uncertainty on the model parameters (see 
next section). 



2.1 Bayesian Evidence 

Although regularisation ensures that the solution given by 
equation (llOfl is well defined, it unfortunately introduces a 
new problem. Regularising a solution reduces the effective 
number of degrees of freedom by an amount that can not be 
satisfactorily determined. Furthermore, applying the same 
regularisation weight to two different models (for example 
different pixellisations of the L—z plane) results in a different 
effective number of degrees of freedom for each model. This 
means the minimum is biased away from the most prob- 
able solution. More crucially, comparison between different 
models cannot be carried out fairly using the statistic. 
Therefore, can not be used to determine the most appro- 
priate pixellisation given the observed data. 

One solution to the problem is to simply not regularise. 
Fortunately, a better solution can be found by turning to 
Bayesian inference and ranking mo d els by their Bayesian 
evidence instead of x^ ■ ISuvu et al. I (|2006f ) derived an ex- 
pression for the Bayesian evidence, e, for the linear inversion 
problem described by equation (|10p . Using the previous no- 
tation, this can be written 



21ne 



X 



In (det[AR]) + In (det[M + AR]) 
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+ Ae^Re + ^ ln(27rCT?) (12) 

i 

with given by equation ((?]). Here, the covariance between 
aU pairs of number count bins has been set to zero (i.e., it 
is assumed that all observed data points are independent of 
ea ch other. For covar iant data, the more general form given 
bv lSuvu et al~ll2006l . would be used). 

The evidence is a probability distribution in the model 
parameters, allowing different models to be ranked fairly to 
find the most probable model. Formally, the evidence should 
be marginalised over all parameters an d the resultin g prob - 
ability used in the ranking. However, ISuvu et al. I (|2006| ) 
noted that a reasonable simphfication is to approximate the 
evidence distribution as a delta function and use the max- 
i mum evidence dir ectly to rank models (see the appendix 
of iDve et al.1l2008l . for more details). We have verified that 
this approximation is still valid for our application and hence 
have adopted it in the present study. 

In our case, the evidence distribution is a function of 
two model parameters; the regularisation weight. A, and a 
parameter, pthrcsh, described in the following section that 
controls the average size of adaptive pixels in the L — z 
plane. Since the evolution depends on the value of these 
two parameters, their uncertainty, as given by the evidence 
distribution, must be included in the overall uncertainty on 
the evolution. For every L — z plane pixel, we compute this 
additional uncertainty, a^, using: 

gI = ^ ej Apthrcsh A(log A) (ej- < e >)^ (13) 

3 

where the sum acts over all points on the grid spanned 
by A and pthrosh at which the evidence, tj, is evaluated. 
Here, Apthresh and A(logA) is the separation between 
grid points along the pthrcsh and log A directions respec- 
tively and < e > is the mean evolution in the pixel. 
Note also that the evidence must be normalised such that 
y^ ^. ej Apthresh A(log A) = 1. This additional error is there- 
fore effectively the scatter in the evolution, weighted by the 
evidence. We compute the overall error in the evolution as 
the quadrature sum of this error and the uncertainty given 
by the appropriate term in the covariance matrix given in 
equation (jlip . 

2.2 Adaptive gridding of the L — z plane 

As we alluded to earlier, the L — z plane can be pixellised 
in a completely general way. This has the advantage that 
smaller pixels can be placed in regions where there are su- 
perior observational data, i.e., higher signal-to-noise (S/N) 
and/or a higher density of data points. The effect of this is 
to better resolve the evolution function whilst maintaining 
an approximately constant level of constraints per pixel over 
the L — z plane. Note that although the method prescribed 
by Peacock & Gull can adjust to the data by limiting the 
series expansion of the evolution function, this alters the 
global resolution over the plane, not in specific regions that 
are better constrained. 

There are several criteria that could be used to control 
the size of pixels in the L — z plane. One example would be 
to make direct use of the density of data points weighted 
by their S/N. Our criterion of choice is to use the covari- 



ance between pixel pairs provided by equation (|lip . Joining 
two pixels together increases the statistical independence 
and hence lowers the covariance of the resulting larger pixel 
with its neighbours. Defining a covariance threshold, pthrosh, 
therefore introduces a criterion whereby pixel pairs whose 
covariance is larger than pthresh are joined together. 

The procedure we have used for adaptively gridding the 
L ~ z plane is as follows: 

1) The L — z plane is uniformly pixellised with a regular 
grid of small rectangular pixels. We tested a range of initial 
grid sizes and found that a 20 x 20 grid is a good compromise 
between having sufficient resolution to properly adapt to the 
varying range of covariances over the L — z plane and being 
of low enough resolution to maintain a high execution speed. 

2) For a given regularisation weight, the covariance matrix 
given by equation (fTTj) is computed. 

3) The lower triangle of the covariance matrix is scanned 
for elements with an absolute value larger than the thresh- 
old covariance. Working from largest to smallest, pixel pairs 
whose covariance exceeds the threshold are joined. If, in 
working down the list, a pixel pair is encountered where 
one of the pixels has already been joined (because it had 
a higher covariance with another pixel), that pixel pair is 
skipped. 

4) New versions of the matrices M and R are computed 
taking into account the new pixellisation. 

Steps 2) to 4) are repeated until the absolute value of all 
pixel covariances lie beneath the threshold pthrcsh. A more 
formally correct procedure would be to recompute the co- 
variance matrix every time a pair of pixels are joined. This 
is considerably slower to execute in practice and we have 
found that the resulting adapted grid does not differ signifi- 
cantly from that produced by the faster procedure outlined 
above. 

2.3 Regularisation 

The regularisation matrix introduced in equation (|10|) arises 
as a result of adding a term to that becomes numerically 
smaller for smoother solutions. In this sense, regularisation 
behaves as a prior, acting to penalise noisy evolution func- 
tions. This is necessary to prevent ill-conditioning in equa- 
tion (|9]). A potential disadvantage is that the evolution may 
not be smooth in reality so that regularisation biases the 
solution. However, the regularisation weight is set by the 
evidence which is maximised at the optimal value of the 
regularisation weight A. If more data points are available, or 
if the data have a higher S/N, the evidence is maximised at a 
lower value of A, allowing greater resolution of any sharp fea- 
tures in the reconstructed evolution. We have investigated 
this to an extent with a synthetic dataset produced using 
an evolution function with a cutoff (see Section |3}. 

The regularisation term added to can be written 

AB = A^ 7?jfce, efe (14) 

where the Rjk are the elements of the regularisation matrix 
R. These depend on the regularisation scheme adopted as 
described below. This means that the solution remains linear 
since the partial derivative of B with respect to the the 
is linear in e. 
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The simplest form of regularisation is zeroth order reg- 
ularisation, where solutions that deviate from zero are pe- 
nalised using Wjk = 5jk so that B = e|. More useful 
forms of regularisation are first order schemes which pe- 
nalise solutions that deviate from a constant and second 
order which penalise solutions that deviate from a gradient. 

For most solutions, as we have found in carrying out 
the work presented herein, first and higher order regularisa- 
tions give very similar results. Whilst we could implement 
first order regularisation, higher orders are generally more 
effective in terms of allowing more freedom in the solution. 
However, conventional higher order schemes are ill-defined 
on our adaptively gridded L — z plane. For example, with a 
second order gradient scheme, pixel triplets required to form 
a finite difference second order derivative would be neither 
co-linear nor parallel. With this in mind, in this work, we 
have opted for a hybrid scheme that uses the difference be- 
tween the evolution in a given pixel j and the sum of evolu- 
tions in all neighbouring pixels k weighted by 

hjk = iQk/aj)Nexp{-y%/r^) . (15) 

Here, Q is the area of the pixel in the L ~ z plane, yjk 
is the separation of the centres of pixels j and k and 
is a normalisation constant set such that hjk ~ 1 

and A'' = 1 when j — k. We set r = 1 in units of pixels. 
The matrix h composed of elements given by equation (jlSp 
relates to the regularisation matrix R via R — h^h. By 
construction, R is non-singular so that the evidence given 
by equation (|12p is always calculable. 

As a final note, all regularisation weights quoted in this 
paper are scaled by the ratio of traces Tr[M]/Tr[R]. This 
normalisation factor means that the regularisation term 
Ae^Re and in equation (|12p are weighted approximately 
equally when the scaled regularisation weight is unity. 

2.4 Evidence maximisation procedure 

There are two non-linear variables which must be ad- 
justed when searching for the maximum evidence. These are 
the regularisation weight, A and the covariance threshold, 
Pthrosh. We also experimented with allowing the initial grid 
resolution to vary as a free parameter but found a strong 
degeneracy with pthresh (if smaller pixels are used, more are 
joined to give a similar final adaptive grid). 

With every trial set of A and pthrosh, the procedure we 
have adhered to is: 

1) Form an adaptive grid (see Section [2. 2 [) . 

2) Form the regularisation matrix R for the current adap- 
tive grid. 

3) Compute the matrix M and vector d as described in 
Section (2] From these, solve for the evolution as given by 
equation pop . We use an LU factorisation routine with it- 
erative refinement. 

4) Compute the covariance matrix as given by equation 
(|lip . given by equation Q and then Ine from equation 

m- 

To find the maximum evidence, we performed a grid 
search over regularly stepped values of logA and pthrosh. The 
nature of the adaptive gridding process is such that the ev- 
idence surface is not perfectly smooth. Despite smoothly 
varying A and pthrcsh, discontinuous jumps in the evidence 



occur when pixels are joined together. By choosing a smaller 
pixel grid, the size of the discontinuities are reduced. Our 
findings indicate that, provided the numerical evaluation of 
the integrals in equations ((3]) and ([6]) is precise, the evidence 
surface is sufficiently smooth with our chosen 20 x 20 grid 
for easy identification the global maximum. In light of this, 
a more efficient automated search for the maximum such 
as a downhill simplex method could be used. In the results 
presented in the next section, we show contour plots of the 
evidence as a function of logA and pthrcsh. 



3 APPLICATION TO SYNTHETIC DATASETS 

In this section, we demonstrate the method by applying it 
to synthetic datasets. These datasets are created from an 
input evolution function which can be compared with the 
reconstructed evolution functions. 

3.1 Synthetic dataset construction 

Three ingredients are required to synthesize a dataset: 

1) Dataset characteristics: In the case of number counts, 
these are the flux limits, in the case of number counts binned 
by redshift, the flux and redshift limits and for background 
flux measurements, the frequencies at which the flux is mea- 
sured. In addition, all three types of dataset must stipulate 
a survey area for determination of uncertainties. 

2) The local luminosity function: To keep our demonstra- 
tion reasonably realistic, we have used a luminosity func- 
tion derived from local galaxies (selected with velocities 
< 30, OOOkm/s) in the Infra-Red Astronomical Satellit e 
(IRAS) Point Source Catalogue of ISaunders et al. I l|2000l ). 
In order to calculate the luminosity of each galaxy at each 
frequency of interest and to compute the bolometric lumi- 
nosity, we estimated an SED. Each SED was fit to the 60 and 
100 nm IRAS fiux and an estimate of the 850 /im fiux pro- 
vided by the tigh t correlat i on between 60, 1 00 and 850 /im 
fiux reported by IVlahakid jPunne fc Ealed ). We refer the 
reader to a forthcoming study where we apply the method 
to real data (Dye et al., in preparation) for more details. 

3) The evolution function: We used two different evolu- 
tion functions. The first is motivated by the findings of 
I Heavens et al. I l|2004l ) and is a smooth function that mono- 
tonically increases with increasing redshift and luminosity. 
We used a bi-linear interpolation in log(l -I- z),logL space 
according to E{L,z) = 1 + xyEmax where x = logj^(j(l -|- 

2)/loglo(l + Zmax), y = (logy,!/) /logjQ (L,„ax /I/min ) and 

Emax = 1000 is the maximum evolution at the point 
{Lmax, Zmax) iu the L — z plauc. The second model evolu- 
tion function is the product of the first and an exponential 
cutoff 1.66 e^~^ which applies at 2 ^ 2. The factor 1.66 en- 
sures that the maximum evolution at 2 = 2 is equal to 1000. 
The monotonic evolution function is defined over ^ 2 ^ 3 
whereas the cutoff function extends to ^ 2; ^ 5. Both 
are defined over 10^° L/Lq s: 10^^. Figure [T] plots the 
functions. 

In the case of number counts and number counts binned 
by redshift, we randomised the synthetic data according to 
Poisson statistics. Datasets were generated assuming either 
0.1 deg^ coverage (the 'low S/N' case) or 10 deg^ coverage 
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Figure 1. The two input evolution functions. Top: The mono- 
tonic function which increases with increasing redshift and lumi- 
nosity. Bottom: The cutoff function formed by multiplying the 
monotonic function by an exponential cutoff cx: e^~^ for z > 2. 
Contours in both plots are in intervals of 100. 



(the 'high S/N' case). Since background flux estimates are 
measured via a difi'erent means in practice, we have random- 
ized the multi-frequency set of synthetic background fluxes 
with a 10% Gaussian error. 

For all four combinations of S/N and evolution function, 
we have generated four different datasets: 

A) Number counts at 850 ^m between the limits 2 mjy (the 
confusion limit at 850/im with a 15 m dish) and 200 mJy in 
six flux bins equally spaced in log flux. 

B) Number counts at 850 /im between the limits 0.1 and 
200 mJy in 10 flux bins equally spaced in log flux. 

C) Number counts at 850 /im between the limits 0.1 and 
200 mJy in 10 intervals binned by 10 redshift intervals 
equally spaced in log(l -I- z). 

D) A set of 10 background fluxes evaluated over the fre- 
quency range 2.4 THz - 300 GHz (125 /im - 1mm) at regu- 
larly spaced intervals in log ( frequency ) . 

Given currently available observing facilities, dataset 
C is somewhat unrealistic. Redshift measurements down to 
fluxes of ~ O.l/xmJy at 850 /.im are well below the confu- 
sion limit of 15m single dish observations and higher resolu- 
tion interferometry presently lacks the sensitivity to gener- 



ate large source samples. However, our choice to include this 
dataset is motivated by the fact that the dataset provides 
an indication of the best possible performance that can be 
expected from the method, thereby setting a benchmark for 
future surveys to aspire to. 

The next section describes application of the recon- 
struction method to different combinations of these datasets. 
In addition to this. Section |3]2]4] outlines further tests of the 
method with a more realistic set of data based on what 
forthcoming Herschel and SCUBA2 surveys are expected to 
deliver. 



3.2 Evolution reconstructions 

We have applied our method to various dataset combina- 
tions to test how its performance depends on data availabil- 
ity and quality. The flrst and second groups of reconstruc- 
tions in Sections I3.2.1l and l3.2.2l applv to datasets generated 
using the monotonic and cut-off evolution function respec- 
tively (see Figure [l}. In the third group presented in Section 
13.2.31 we compare recovery of the evolution with high S/N 
versions of datasets in the first and second groups. Finally, 
the last set of reconstructions in Section 13.2.41 are based on 
synthetic Herschel and SCUBA2 datasets. 

3.2.1 Monotonic evolution 

The first group of reconstructions are shown in Figure [S] 
These apply to the four dataset combinations A, A+D, B 
and C (see Section [3. If) generated using the monotonic evo- 
lution function in the low S/N case. The top row in the figure 
plots the coverage of each dataset in the L — z plane. For 
each of the four test reconstructions shown, we have binned 
the input evolution function to the adaptive grid determined 
in each case to facilitate comparison with the recovered evo- 
lution. The binned input and recovered evolution is plotted 
in the second and third rows of the figure respectively. Com- 
paring the two clearly shows that recovery of the input evo- 
lution becomes increasingly more faithful as the coverage of 
the datasets in the L — z plane increases. The average size 
of adaptive pixels selected by the maximum evidence de- 
creases with increasing coverage as a result of the increased 
constraining power of the data. Similarly, as the evidence 
contour plots at the bottom of the figure show, less regular- 
isation of the evolution function is necessary when the data 
coverage is more extensive. Notice also that the covariance 
threshold at which the evidence is maximised is lower when 
the constraints provided by the data are greater. Finally, an 
important point to note is that in all cases, the significance 
of the residuals is distributed approximately normally with 
unit standard deviation. This means that the uncertainty on 
the evolution derived by our procedure is an accurate repre- 
sentation of where the recovered evolution truly differs from 
the input evolution. 

The evidence contours show that a mild degeneracy ex- 
ists between the covariance threshold and the regularisation 
weight. The degeneracy arises because a higher regularisa- 
tion weight imposes stronger constraints on pixels. This has 
the effect of reducing their average covariance as quantified 
by equation (|lip . Therefore, at a fixed covariance threshold, 
as the regularisation weight increases, the average pixel size 
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Figure 2. Application of the reconstruction method to four different combinations of low S/N datasets generated using the monotonic 
evolution function. Reading from left to right, columns correspond to datasets A, A+D, B, C (see Section 13.11 1. Reading from top 
to bottom, each row shows 1) coverage of input datasct(s), 2) input evolution function binned to the adaptive grid determined in the 
reconstruction (contours show the smooth underlying function plotted in Figure^, 3) reconstructed evolution function, 4) the significance 
of the reconstructed evolution function, 5) the significance of the absolute residuals between the binned input and reconstructed evolution 
functions, 6) evidence contours in the plane spanned by covariance threshold and regularisation weight. Contours are stepped in four 
intervals of A Ine = 1 away from the maximum. 
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decreases. To maintain an optimal adaptive grid (i.e., an 
optimal average pixel size) at higher regularisation weights, 
a lower covariance threshold must be used and vice versa. 
In fact, the inclination of the evidence contours is parallel 
to lines of constant numbers of pixels. The maximum ev- 
idence therefore selects a specific subset of adaptive grids 
with approximately the same number of pixels from the set 
of all possible grids. The contours do not extend indefinitely 
because at the low covariance threshold end, the higher reg- 
ularisation reduces the dynamic range of covariances over 
the L — z plane. This results in an adaptive grid with more 
similarly sized pixels, i.e., the grid is not able to optimally 
adapt to the variation in data coverage. The reverse is true 
at the high covariance threshold end, where the lower reg- 
ularisation results in a larger range of covariances to the 
extent that some pixels become too large while some remain 
too small to properly adapt to the data coverage. 

Considering each reconstruction in turn in Figure [21 it 
is apparent that number count data alone are unable to re- 
cover the input evolution to a satisfactory accuracy. In the 
case of dataset A shown in the first column (with number 
counts that reach down to 2 mjy and no background flux 
measurements), the evolution is biased low at high L, high 
z and biased high at high L, low z and at low L, high z. 
This biasing occurs because the evolution function has been 
'smeared' over the L — z plane by two effects. The flrst effect 
causes the evolution to be smeared along the lengths of flux 
bins. This is a direct result of the fact that the number of 
observed galaxies in a given bin can be reproduced by many 
different evolution proflles along the bin's length, a point we 
will return to below. The second effect is that the regularisa- 
tion effectively smooths the evolution. Regions of the L — z 
plane where there are no constraints from the data are only 
constrained by the regularisation. In these regions, the evo- 
lution is extrapolated from neighbouring data-constrained 
regions, by an amount governed by the regularisation weight. 
This is a useful feature which we alluded to in the introduc- 
tion: The ability to extrapolate is one of the key advantages 
of model-based approaches to measuring evolution. 

It is apparent that the residuals in the flrst column of 
Figure [2] show that the most significant deviation from the 
input evolution occurs not where data constraints are lack- 
ing, but at high L, high z where the evolution is strongest. 
This is again due to the smearing which causes a larger abso- 
lute difference when the evolution is stronger. In this sense, 
it could be argued that the method has performed more re- 
liably in regions of the L — z plane not covered by data. Of 
course this is primarily because the method assigns a much 
higher uncertainty to areas of zero coverage. The ratio of 
input to recovered evolution in areas of zero data coverage 
shows a stronger deviation from unity than areas that are 
covered by data. 

Despite the biasing, the number count data still suc- 
cessfully detect a significant increase in the evolution from 
low L, low z to high L, high z. At first sight, this seems 
incredible given that number count data contain nothing 
more than the number of galaxies between two fiux limits. 
In other words, a single flux bin gives only a measure of the 
average evolution within it and nothing about the variation 
of the evolution. The same measured number of galaxies be- 
tween two fiux limits can be reproduced with a wide range 
of different evolution functions. Although true of an indi- 
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Figure 3. Synthetic 850 ^tm cumulative number counts generated 
with five different evolution functions. Four of these are the four 
90-degree rotations of the monotonic evolution function plotted 
in Figure [T] such that the maximum lies at each of the four cor- 
ners of the L — z plane. The fifth evolution function is described 
by a Gaussian with a fuU-width-at-half-max of A log L = 1 and 
l\z = 1, reaching the same maximum evolution but peaking in 
the centre of the L — z plane. The counts arising from this fifth 
evolution function are shown by the broader grey line. 



vidual bin, this degeneracy becomes greatly reduced when 
other fiux bins are added to the constraining data, especially 
when the solution is regularised. 

FigureOillustrates the constraining power of pure num- 
ber count data. We have plotted synthetic 850 /xm number 
counts generated using five different evolution functions cov- 
ering 10^" ^ L/Lq ^ 10" and < z ^ 3 as in Figure [21 
Four of these are monotonic and are simply the four 90- 
degree rotations of the previous monotonic function so that 
the maximum evolution is located at each of the four corners 
of the L — z plane. The fifth evolution function was chosen to 
be a Gaussian with a fuU-width-at-half-max of A log L = 1 
and Az — 1, reaching the same maximum evolution but 
peaking in the centre of the L — z plane. 

The figure shows that the four different monotonic 
evolution functions give rise to four unique number count 
datasets. All four datasets remain unique if any one under- 
goes the renormalisation n(/) — > OLn{f) where a is a scale 
factor. This means that none of the four evolution functions 
can be renormalised to create a set of number counts that 
match those obtained using any of the other three evolution 
functions. Reversing this argument, the number count data 
can therefore distinguish between the four evolution func- 
tions. However, the Gaussian evolution function produces 
number counts which are very similar under renormalisa- 
tion to the monotonic function that peaks at high L, high 
z (compare the broader grey line with the thin continuous 
black line in Figure [Sj. This occurs because translating the 
Gaussian function by Az = -|-1.5 along the direction of the 
850 pim fiux bins (i.e., so that the peak lies at z = 3), results 
in an evolution function that is very similar to the high L, 
high z monotonic function. This serves to demonstrate the 
insensitivity of flux bins to variations in evolution along their 
length. The four rotated monotonic functions can be distin- 
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guished because they have moments perpendicular to hues 
of constant flux in the L~ z plane. Indeed, the reconstructed 
evolution obtained from the number counts generated using 
the Gaussian evolution function is very similar to the mono- 
tonic function which peaks at high L, high z, apart from a 
normalisation offset. Number counts can therefore not di- 
rectly measure evolution cut-offs in redshift. We return to 
this point below. 

The second column in Figure [5] shows the method ap- 
plied to the combination of datasets A and D (i.e., 850 /xm 
number counts down to 2mJy and 10 measurements of the 
background). Compared with just the number counts shown 
in the first column, it is clear that the reconstruction is more 
accurate. The added constraints from the background flux 
reduce the level of smearing. This is particularly apparent at 
low luminosity and high redshift where the biasing is much 
smaller. The additional constraints also enable a higher res- 
olution reconstruction at low L, low z and result in a lower 
required level of regularisation. 

In the third column of Figure [2] we show how well the 
evolution is recovered with dataset B which comprises pure 
number count data again, but down to a lower flux limit of 
0.1 mjy. In this case, the flux bins cover almost the entire L— 
z plane. The figure shows that a similar, if not slightly higher 
level of biasing in the recovered evolution occurs, compared 
to that seen in the second column, due to the same smearing 
effects. However, there is an obvious improvement compared 
to the shallower number counts shown in the first column. 
This is apparent at low L where the shallow counts permit 
a lower resolution reconstruction due to lack of coverage. 
Notice also that the solution requires a much lower level of 
regularisation. 

A point to note is that the reconstructed evidence in 
regions of the L — z plane constrained by data are not no- 
ticeably biased by regions that are not constrained by data. 
For instance, comparing the first and third columns of Fig- 
ure [2] in the region of the L — z plane covered by the bright- 
est six flux bins shows that the recovered evolution and the 
significance of the residuals are very similar. Therefore, in- 
terpolation into areas lacking data coverage can be relied 
upon without detriment to the overall solution. 

The fourth column of Figure [2] corresponds to the case 
where redshift information is included. Not surprisingly, the 
addition of redshifts improves the reconstruction dramati- 
cally by providing direct measurements of the evolution at 
specific points in the L — z plane. The dataset allows re- 
construction of the evolution on a much higher resolution 
grid. The larger pixels that still remain at high luminosities 
are due to increased Poisson noise caused by a significantly 
lower number density of galaxies (a common feature of most 
of the results shown in this paper). Despite the increased 
resolution, the significance of the reconstructed evolution 
is 3 times higher compared to the cases where redshifts 
are not present. The solution is again biased most at high 
L, high z where the evolution peaks and smearing has the 
largest effect, but this is less extensive and much reduced in 
magnitude. 

3.2.2 Evolution with a cut-off 

The second group of reconstructions apply to the four 
dataset combinations A, A-fD, B and C (see Section IS.lf) 



generated in this case using the cut-off evolution function 
with low S/N (i.e. assuming an area of 0.1 deg^). The top 
row in Figure |3] shows the coverage of the datasets over the 
L — z plane which this time extends further in redshift out 
to 2 = 5. 

In terms of the variation of characteristics such as res- 
olution, residuals, significance of the reconstructed evolu- 
tion, optimal regularisation weights and optimal covariance 
thresholds, the reconstructions exhibit very similar traits to 
those seen with the previous group of datasets generated us- 
ing the monotonic evolution function. However, there is one 
obvious and significant difference. Although the monotonic 
evolution function appears to be fairly well reconstructed 
without redshift information, the cut-off evolution function 
can not be properly recovered unless redshifts are included. 
Smearing along the length of flux bins removes any sign of 
the cut-off, pushing the strongest evolution to the highest 
redshift. This is the exact same behaviour observed previ- 
ously with the ad-hoc Gaussian evolution function. We thus 
iterate again that redshift data is essential to measure any 
breaks or discontinuous features in the evolution. 

We have investigated where, in relation to the break, 
redshift data most effectively constrains the break. Our find- 
ings show that redshift measurements lower than the break 
make little or no difference to identifying the location or 
even the existence of a break. Only redshift measurements 
at the break and beyond are able to constrain its location. 



3.2.3 High S/N reconstructions 

We have repeated some of the reconstructions carried out 
previously but with high S/N datasets (i.e. where an area of 
10 deg^ has been assumed) . We generated the two extremes 
in datasets considered, namely the shallow number count 
only dataset and the deeper number counts binned by red- 
shift, using both the monotonic and the cut-off evolution 
function. Figure [5] shows the results. 

Comparing with the corresponding low S /N reconstruc- 
tions, the most striking difference is the improved resolution 
and higher significance of the reconstructed evolution. The 
rib-like structures present in the evolution significance (also 
seen to a lesser degree in the low S/N case) occur because 
of the varying number of constraints per pixel which range 
from one to four data bins depending on their alignment 
with the L — z plane pixel grid. 

An increase in the accuracy (i.e., a reduction in the av- 
erage absolute size of the residuals) is also clear. Although 
the range of residual significances is approximately the same 
between high and low S/N datasets, the extent of the most 
significant discrepancies is less in the high S/N cases. In fact, 
the distribution of residual significances is slightly narrower 
than a normal distribution with unit Icr width, implying 
that the errors are probably slightly over-estimated. If the 
residuals were randomly scattered over the L — z plane, then 
this would imply that the evolution had been recovered to 
the best accuracy possible. However, the fact that the resid- 
uals form coherent patterns resembling the evolution func- 
tion indicates that the reconstruction is not perfect. There 
are still low- level biases present, an inescapable consequence 
of regularisation. 
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Figure 4. As Figure [2l but applied to low S/N datasets generated using the cut-ofF evolution function. 
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Figure 5. Application of the reconstruction method to high S/N datasets. The first and second columns correspond respectively to 
datasets A and C generated with the monotonic evolution function and the third and fourth columns to datasets A and C generated 
with the cut-off evolution. Row descriptions arc the same as those in Figure [J] 
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3.2.4 Forthcoming survey predictions 

The reconstructions previously considered give a good indi- 
cation of how well different dataset types are able to con- 
strain the evolution. From these simulations, it is clear that 
redshift data is vital to properly recover the input evolution. 
However, as we discussed in Section 13.11 the extent of the 
redshift data coverage in the L — z plane was overly opti- 
mistic. Furthermore, in practice, observational data gath- 
ered for the purpose of measuring evolution will most likely 
be heterogeneous, comprising redshift and number count 
measurements from different studies observed at different 
wavelengths and to different depths. 

In this section, we investigate the performance of the 
method with synthetic data generated to mimic a more 
realistic collection of data, typical of what might be ex- 
pected from forthcoming submm surveys. Two of the largest 
datasets will be delivered by Herschel and SCUBA2. In the 
case of Herschel, we synthesized data based on the Spec- 
tral and Photometric Imaging Receiver (SPIRE) which op- 
erates in three passbands with central wavelength 250, 350 
and 500 /im. For SCUBA2, we synthesized data at its longer 
wavelength channel at 850 ^m. 

In terms of redshifts, we have assumed that measure- 
ments can be made down to the flux confusion limit at 
eac h wavelength. For SPIRE, th ese limits, as estimated 
by iLagache. Dole fc Pugetl (|2003l ). are approximately 10, 
20 and 20mJy at 250, 350 and 500 /im respectively. The 
confusion limit for SCUBA2 at 850 ^m, is ~ 2mJy (e.g., 
iHughes et al. |[l998l ). In terms of number counts, we have 
assumed that a flux limit 20 x lower than the confusion limit 
can be achieved using technique s such as the so-called 'P (D)' 
statistical method (see, for e.g.. IPatanchon et al. |[2009l . and 
references therein) or to a lesser extent, exploitation of grav- 
itational lens magnification (e.g.. Small, Ivison & Blain 1997 
or more recently, Knudsen et al. 2006). 

We tested the method with the synthetic data gener- 
ated using the two different evolution functions as before. 
The top four panels of Figures |6] and [7] show the coverage of 
the datasets in the case of the monotonic and cut-off evolu- 
tion function respectively. The redshift data cover the flux 
ranges 10 - 1000 mjy, 20 - 1000 mjy, 20 - 500 mjy and 2 
- 200 mJy at 250, 350, 500 and 850 /im respectively in six 
equal intervals in log(L) and 10 equal intervals in log(l -I- 2). 
Similarly, the number count data cover four flux bins log- 
arithmically spaced between the ranges 0.5 - lOmJy, 1 - 
20mJy, 1 - 20mJy and 0.1 - 2mJy at 250, 350, 500 and 
850 fj,m respectively. We assumed an area of 1 deg^ . 

The bottom four panels in Figure[6]show the reconstruc- 
tion results for the monotonically evolving data. Despite a 
10 times reduction in area and lower redshift data coverage 
of the L — z plane compared with the high S/N case (shown 
in the second column of Figure [SJ, the input evolution is 
recovered with a comparable accuracy. At high luminosities, 
the average significance of the residuals is very similar to 
that in the high S/N case. This is not too surprising since 
the relative reduction in area is offset by the more numerous 



^ We have not considered the two shorter wavelength passbands 
of the Photodetector Array Camera and Spectrometer on board 
Herschel since the corresponding flux bins only cover a small re- 
gion of our L — z plane. 
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Figure 6. Application of the reconstruction method to sim- 
ulated HersciieJ/SPIRE (250, 350 and 500 /xm) and SCUBA2 
(850 fim) data generated using the monotonic evolution func- 
tion over 1 deg^ . The top four panels show the coverage of each 
dataset. In each case, redshift data extends down to the confu- 
sion limit and number counts reach down another factor of 20 in 
flux (see main text). The bottom four panels show the binned in- 
put evolution, the reconstructed evolution, the signiflcance of the 
reconstructed evolution and the absolute significance of the resid- 
uals. The evidence in this case is maximised at Pthresh = 2800, 
logio A = -2.2. 



redshift data across the different wavelengths. Although the 
residuals fare worse at low luminosity where redshift data 
are completely lacking, they are smaller than the residuals 
that arise when no redshift data are present anywhere on 
the L — z plane. For example, the residuals at low L, high 
z in the third column of Figure [2] are higher and evaluated 
over larger pixels. The difference with the reconstruction 
using multi-wavelength data is that the variation of lumi- 
nosity with redshift for a given flux changes as a function 
of observed wavelength. This is a direct result of the con- 
version between monochromatic and bolometric luminosity 
which depends on the SED assumed. The effect is that lines 
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Figure 7. Application of the reconstruction method to simulated 
HerscheJ/SPIRE (250, 350 and 500 ^m) and SCUBA2 (850 iira) 
data generated using the cut-off evolution function over 1 deg^ . 
Panels are as described in Figure [B] The evidence in this case is 
maximised at Pthrcsh = 2700, logio A = —2.8. 



of constant flux, and therefore flux bins, in the L — z plane 
corresponding to two different observed wavelengths can in- 
tersect each other. In the case of the reconstruction shown 
in Figure [H some of the flux bins without redshift informa- 
tion cross other bins corresponding to different wavelengths 
where redshifts are present. Since data binned by both flux 
and redshift give a direct measure of the evolution at a spe- 
ciflc location in the L — z plane, this provides additional 
constraints on how the evolution must vary along the length 
of a flux-only bin. 

Figure[7]shows how the method performs in reconstruct- 
ing the simulated iIersciieJ/SCUBA2 data generated using 
the cut-off evolution function. As in the case of the mono- 
tonically evolving data, the residuals exhibit a similar level 
of variation. However, unlike the monotonic case, the res- 
olution of the recovered evolution is slightly lower than in 
the high S/N case. This is especially true at the highest 
redshifts where the evolution is considerably weaker than in 



the monotonic case. This results in much smaller numbers 
of galaxies per unit area of the L — z plane and hence in- 
creased Poisson noise and covariance. The adaptive gridding 
procedure therefore increases the pixel size in these regions. 



4 SUMMARY 

In this paper, we have presented a new method for measuring 
evolution of the galaxy luminosity function. The method 
computes the evolution as a discretised function of redshift 
and bolometric luminosity. The advantage brought about 
by this discretisation is that the evolution can be solved 
linearly, ensuring the best fit solution is always found with 
a single, efficient matrix inversion. Formulating this as a 
linear problem also brings the additional advantage that the 
uncertainty in the reconstructed evolution can be obtained 
with a few simple extra computational steps. Furthermore, 
by introducing linear regularisation, the evolution can be 
extrapolated into regions of the L — z plane where data are 
lacking. 

We have also developed a procedure for adaptively pix- 
ellising the evolution function. This allows the evolution to 
be reconstructed with higher resolution in regions of the 
L — z plane where data constraints are stronger and lower 
resolution where constraints are weaker. The procedure uses 
the data to automatically set the optimal pixellisation via 
maximisation of the Bayesian evidence. 

We have applied the method to a range of synthetic 
datasets constructed with two different input evolution func- 
tions; one rising monotonically with redshift and luminosity 
and the other incorporating a redshift cutoff. Comparing the 
reconstructed evolution with the input evolution provides a 
means of testing how well the method performs with varying 
dataset type, coverage and signal-to-noise. Our findings in- 
dicate that redshift measurements are essential to locate the 
presence of a cutoff in the evolution and that these must lie 
beyond the cutoff. Number count data alone allow for a sur- 
prisingly reasonable reconstruction if the evolution function 
is known to be monotonic, but falsely indicate monotonic 
behaviour when a redshift cutoff is present. 

We have made predictions of the degree to which forth- 
coming submm surveys carried out with new instruments 
(SCUBA2 and Herschel) will allow evolution in the submm 
GLF to be determined. These simulations show that com- 
bining a mixture of number count data from both facilities 
and including measurements of redshifts down to the source 
confusion limit will allow for a very well resolved and reli- 
able measurement of evolution. In particular, the predictions 
have demonstrated that much improved constraints are pro- 
vided by pure number count data measured at a combina- 
tion of different wavelengths, compared to counts measured 
at just one wavelength. This is a result of the cross-linking 
of flux bins on the L — z plane due to the conversion be- 
tween monochromatic and bolometric luminosity which is a 
function of both redshift and wavelength. 

In terms of applying the method to real data, there 
are a few additional complexities that must be considered. 
One complexity is that real data are typically incomplete. 
This will result in the reconstructed evolution being a lower 
limit on the real evolution. Indeed, our simulations indi- 
cate that, due to regularisation, the reconstructed evolution 
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tends to be underestimated even when the data are com- 
plete, although by an amount which is within the derived 
uncertainties. However, depending on the magnitude of the 
incompleteness, this may give rise to an underestimate that 
is not within the derived evolution error budget. In addi- 
tion, depending on the level of regularisation, severe incom- 
pleteness in one dataset could significantly affect the recon- 
structed evolution in other areas of the L — z plane covered 
by near-complete data. 

Another consideration that must be made with real data 
is regarding the SED used for converting between monochro- 
matic and bolometric luminosity. Clearly, if this is a poor 
match to the SEDs of the galaxies in a given dataset, then 
the assumed location of that dataset in the L — z plane 
will be inaccurate and give rise to a systematic error in the 
recovered evolution. In the submm, this will be larger for 
fluxes near the peak of the observer- frame SED, but on av- 
erage, should be a relatively small effect. Of course, instead 
of using a single SED as we have done in our simulations, 
source-specific SEDs, or dataset-specific SEDs could be used 
in practice. 

In this paper, we have demonstrated the behaviour of 
the reconstruction method with a combination of different 
datasets and evolution functions. Although we have con- 
centrated on key characteristics, there are many more that 
could have been tested. Since an exhaustive investigation 
of all configurations that might be encountered in practice 
is well beyond the scope of this work, a sensible approach 
would be to carry out simulations when applying the method 
to datasets that differ greatly from those tested here. In this 
way, any potential systematics that might arise from the 
data or any non-uniqueness in the reconstructed evolution 
could be quantified. 

There are also several ways in which the method 
could be developed. For example, the adaptive pixellisation 
scheme is relatively simple and could be enhanced to offer 
greater adaptability. Another possible enhancement might 
be improving the use of redshift information. In this work, 
we incorporated redshifts by crudely binning number counts 
by redshift. This lowers the constraining power of the red- 
shift data. With this in mind, we have begun development 
of a more sophisticated scheme that includes redshift infor- 
mation on a more efficient source by source basis. This will 
be presented in forthcoming work. 

We are about to enter a new era in submm cosmology 
with the arrival of new instruments such as Herschel and 
SCUBA2. Surveys conducted with these instruments will 
give an increase in the number of detected submm sources 
of several orders of magnitude compared to existing surveys. 
Their significantly improved sensitivities will allow much 
wider ranges in luminosity and redshift to be explored. Com- 
bined with new methods to measure evolution, over the next 
few years, this will bring about a revolution in our overall 
understanding of how galaxies form and evolve. 

Acknowledgements 

SD acknowledges support by the Science and Technolo- 
gies Facilities Council. We would like to thank the reviewer 
of this paper. Professor Steve Phillipps, for his helpful com- 
ments and suggestions. 



REFERENCES 

Dunlop, J. S. & Peacock, J. A., 1990, MNRAS, 247, 19 

Dye, S., 2008, MNRAS, 389, 1293 

Eales, S. A., et al., 2009, ApJ, 707, 1779 

Fixsen, D.J., Dwek, E., Mather, J.C., Bennet, C.L., Shafer, 

R.A., 1998, ApJ, 508, 123 
Heavens, A. F., Panter, B., Jimenez, R., Dunlop, J., 2004, 

Nature, 428, 625 
Hughes, D. H., et al., 1998, Nature, 394, 241 
Knudsen, K. K., Barnard, V. E., Van der Werf, P. P., 

Vielva, P., Kneib, J. -P., Blain, A. W., Barreiro, R. B., 

Ivison, R. J., Small, I., Peacock, J. A., 2006, MNRAS, 

368, 487 

Lagache, C, Dole, H. & Puget, J. -L., 2003, MNRAS, 338, 

555 

Patanchon, C, et al., 2009, ApJ, in press. larXiv:0906.0981l 
Peacock, J. A., 1985, MNRAS, 217, 601 
Peacock, J. A. & Gull, S. F., 1981, MNRAS, 196, 611 
Robertson, J. C, 1978, MRAS. 182, 617 
Robertson, J. C, 1980, MNRAS, 190, 143 
Saunders, W., et al., 2000, MNRAS, 317, 55 
Small, I., Ivison, R. J. & Blain, A. W., 1997, MNRAS, 490, 
L5 

Suyu, S. H., MarshaU, P. J., Hobson, M. P., Blandford, R. 

D., 2006, MNRAS, 371, 983 
Vlahakis, C, Dunne, L. & Eales, S. A., 2005, MNRAS, 364, 

1253 

WaU, J. v.. Pope, A. & Scott, D., 2008, MNRAS, 383, 435 
Warren, S. J. & Dye, S., 2003, ApJ, 590, 673 



